%% mf_analyze, print and plot results fo simulation

%% Analyzing output of xsecton model
% Plot /shows the results of the smulation for a given stress period and
% time step to be set by the user
% TO 090223 111212

clear variables
close all

film=1;             % if film==1, make avi file else not

%% load model name and basename contained in name.mat
load name           % Load 
load(basename);     % this yields the stored matrices that make up the grid and the model
load underneath     % Load data from mf_adapt
load setup          % Load setup data from SETUP sheet in basename.xsl

load CurrentSeriesNumber.mat
n=CurNum;           % set series for runall.m

%% Start
[NRow,NCol,NLay]=size(IBOUND); IL=find(IBOUND==-2); IR=find(IBOUND==-4);

%% compute coordinates of axes and cell centers

[xGr,yGr,zGr,xm,ym,zm,Dx,Dy,Dz,Nx,Ny,Nz]=modelsize3(xGr,yGr,zGr);

%% load the unformatted files with the heads, the concentration and the budget terms

H=readDat([basename,'.hds']);
H=maskHC(H,IBOUND & H(1).values<1000);   % read heads and mask inactive cells with NaN

C=readMT3D('MT3D001.UCN');
C=maskHC(C,IBOUND);   % read concentrations and mask inactive with NaN
C=maskHC(C,ICBUND);   % read concentrations and mask inactive with NaN

B=readBud([basename,'.bgt'],'FLOWRIGHTFACE');  B=maskHC(B,IBOUND); % read budgetfile and mask

% for i=1:length(B)
%   for j=1:length(B(i).term)
%       B(i).term{j}(IBOUND==0 | H(1).values>1000)=NaN;
%   end
% end

B=mf_Psi(B);
dPsi=20;
crange=ContourRange(B,dPsi,'','Psi');
minPsi=min(crange);
maxPsi=max(crange);

% Stream function mask, to not plot stream lines above sand body
PSIMASK=NaN(Nz+1,Nx-1);
PSIMASK(inpolygon(ones(size(zGr(:)))*xGr(2:end-1),zGr(:)*ones(size(xGr(2:end-1))),xSand,zSand))=1;

concs= -0.0001:0.05:1.05;           % Conc is always between 0 and 1 
chead=0:0.05:17;                    % Heads is always between base and top of model
cpsi =minPsi:(maxPsi-minPsi)/50:maxPsi; % Psi-range is computed above


%% Photo to plot on and create subplots and setup videofile
foto='IMG_20120131_134330.jpg';     % first photo
load(['.' filesep 'Photos' filesep '2' filesep 'photostruct.mat']); % yields struct d
i=strmatchi(foto,{d.name});

figure('Visible','Off')             % actual figure window is not displayed, faster
fig=figure();                       % create handle
Ax1=subplot(2,1,1);
Ax2=subplot(2,1,2);

ImageModel3(Ax1,['Photos' filesep '2' filesep foto],[d(i).uMdl],[d(i).vMdl],[d(i).DX],[d(i).DY]);
ImageModel3(Ax2,['Photos' filesep '2' filesep foto],[d(i).uMdl],[d(i).vMdl],[d(i).DX],[d(i).DY]);


ImageModel2(Ax1,['Photos' filesep '2' filesep foto],[d(i).uMdl],[d(i).vMdl],[d(i).DX],[d(i).DY]);
ImageModel2(Ax2,['Photos' filesep '2' filesep foto],[d(i).uMdl],[d(i).vMdl],[d(i).DX],[d(i).DY]);


% maxfig(gcf,1);        % To create a high resolution film

dates(1:length(d))=0;   % ugly method to get array of dates
for i=1:length(d)
    dates(i)=d(i).datenum;
end


aviobj=VideoWriter([basename '_cycle_' sprintf('%d',n)]);   % create different video for each serie
aviobj.FrameRate = 10;
aviobj.Quality = 100; %1=Lowest quality, 100=Max quality
open(aviobj);

subplot(2,1,1)
set(gca,'clim',[min(concs) max(concs)]);

%% Freshwater stored
    Quantity=[0;0];
    ConcentrationWell=[0;0];
    ConcentrationMonitorPoint1=[0;0];
    MaxQuantity=dy*dx*SETUP_XWALL(n)*dz*SETUP_ZWALL(n)*POREFF(n)/100;
    MaxAllowedC=125; %[mg/l]
    ExtractionStart=90;  %Time start extracting[min]
    dummy=0;
    dummy2=0;
    MaxUsableCol=0;
    SizeTooSmall=0;
    FISRTCOL=0;


for i=1:length(C);

    A=XS(C(i).values); %rename concentration values
        Quantity=[Quantity, [XS(C(i).time);length(find(A(1:izWall,1:ixWall)<(MaxAllowedC/29000)))*dy*dx*dz*POREFF(n)/100]]; %[cm3]
        FIRSTCOL=A(:,1);
        if (FIRSTCOL>=0.05 & FIRSTCOL<=0.95)
        [mu(i),sig(i)]=normfit(diff(FIRSTCOL(FIRSTCOL>0.05 & FIRSTCOL<0.95)));
        end
        %%  Concentration Monitor Point 1 || Extra concentration monitoring
%         ConcentrationMonitorPoint1=[ConcentrationMonitorPoint1, [XS(C(i).time);A(izWall, ixWall)*29000]]; %Concentration at Well
%     if (ConcentrationMonitorPoint1(end,end)<=0.95*29000 && dummy2==0);
%         dummy2=1;
%         SizeTooSmall=1;
%         SizeTooSmallTime=ConcentrationMonitorPoint1(1,end);
%     end
    %%  Concentration Well
        ConcentrationWell=[ConcentrationWell, [XS(C(i).time);A(izMilk, ixMilk)*29000]]; %Concentration at Well
    if (ConcentrationWell(end,end)>=MaxAllowedC && XS(C(i).time)>=ExtractionStart && dummy==0);
            APrev=XS(C(i-1).values); %rename previous concentration values
            MaxUsableQuantity=[XS(C(i-1).time);length(find(APrev(1:izWall,1:ixWall)<0.0052))*dy*dx*dz];
            CurConcentration=ConcentrationWell(:,end);
            MaxUsableCol=i+1;
            dummy=1;
    end
    %%  Plot info
    if i==1
        axes(Ax1);
        [~,hc]=contourf(xm,squeeze(zm),XS(C(1).values),concs,'g');
%         [~,hh]=contour (xm,squeeze(zm),XS(H(1).values),chead,'r');
        [~,hp]=contour (xGr(2:end-1),squeeze(zGr),B(i).Psi,cpsi,'w');

        xlabel('x [cm]');                 % set x-label
        ylabel('z [cm]');                 % set y-label
        set(gca','color',[1 1 0.8])      % light blue background color

    else
    %%  Find photo
    axes(Ax2);
        CurTime=addtodate(d(1).datenum,C(i).time,'minute');
        ArrayPhoto=find(dates<CurTime);
        if isempty(ArrayPhoto)
            IndexArrayPhoto=1;
        else
            IndexArrayPhoto=ArrayPhoto(end);
        end
        ImageModel2(Ax2,['Photos' filesep '2' filesep d(IndexArrayPhoto).name],[d(IndexArrayPhoto).uMdl],[d(IndexArrayPhoto).vMdl],[d(IndexArrayPhoto).DX],[d(IndexArrayPhoto).DY]);
        daspect([1 1 1]);
        title(sprintf('%s Experimental Setup; picture taken at %s',basename,datestr(d(IndexArrayPhoto).datenum)));
        
        xlabel('x [cm]');                 % set x-label
        ylabel('z [cm]');                 % set y-label        
        
    axes(Ax1);
        period=C(i).period;
        tstp  =C(i).tstp;                   % choose a timestep period present in the file
    
%         set(hh,'zdata',XS(H(i).values));
        set(hc,'zdata',XS(C(i).values));   
        set(hp,'zdata',B(i).Psi);
%         hold on;
        
    end
    axes(Ax1);
        set(get(hc,'children'),'edgecolor','none');
%         set(get(hc,'children'),'edgecolor','none','facealpha',0.15);

        title(sprintf('%s %d th setup, time = %g min\n%s',basename,n,C(i).time,datestr(addtodate(d(1).datenum,C(i).time,'minute'))));

    drawnow;
    if (film==1 && (dummy==0 || REUSENEXTCONC(n)==0)==1);
        writeVideo(aviobj,getframe(fig));
    end        
end
axis tight

if film, close(aviobj); end %save movie.mat F

% Recovery efficiency

% Injection=[Quantity(1,:);0,QMilk'.*(1000/60)*MULTIMILK(n)];  % Injection[cm3/delta t],QMilk[l/h],Diffcon[min]
Injection=[Quantity(1,:);0,QMilk'.*diff(ConcentrationWell(1,:))*MULTIMILK(n)];  % Injection[cm3/delta t],QMilk[cm3/min],Diffcon[min]
% Injection=[Quantity(1,:);0,QMilk'./diff(ConcentrationWell(1,:))*(1000/60)*MULTIMILK(n)];  % Injection[cm3/delta t],QMilk[l/h],Diffcon[min]

[m]=find(cumsum(Injection(2,:))==max(cumsum(Injection(2,:))),1,'first');  %Col where injection stoppes
[o]=find(cumsum(Injection(2,:))==max(cumsum(Injection(2,:))),1,'last');  %Col where dejection startes
QInjected=[cumsum(Injection(2,1:m)), ones(1,length(C)+1-m)*max(cumsum(Injection(2,1:m)))];
QRetreived=[zeros(1,o),cumsum(Injection(2,o:length(C)))];
if (ExtractionStart>=ConcentrationWell(1,o) || m==0 || o==0)
error('Please set ExtractionStart time=%d[min], in mf_analyzepretty on line 84.\n',ConcentrationWell(1,m));
end
 
RE=(QRetreived./QInjected)*-100;
%% Graphs
time=1:ConcentrationWell(1,end);

figure(2);
plot(Quantity(1,:),Quantity(2,:),time,MaxQuantity,Quantity(1,:),QInjected,Quantity(1,:),QRetreived, 'LineWidth',1, 'LineSmoothing','on')
line([0 ConcentrationWell(1,end)],[MaxQuantity MaxQuantity],'Color','black');
title('Amount of freshwater stored in Zoetwatervat, with concentration < 150[mg/l]');
xlabel('Time[min]');
ylabel('Quantity[cm3]');
legend('Quantity in zoetwatervat','Absolute max. quantity zoetwatervat','Q Injected cummulative','Q Retreived cummulative')

figure(3);
plot(ConcentrationWell(1,:),ConcentrationWell(2,:),ConcentrationMonitorPoint1(1,:),ConcentrationMonitorPoint1(2,:),time,MaxAllowedC, 'LineWidth',1, 'LineSmoothing','on')
line([0 ConcentrationWell(1,end)],[MaxAllowedC MaxAllowedC],'Color','black');
title('Concentration at intake');
xlabel('Time[min]');
ylabel('concentration [mg/l]');
legend('Concentration at well','Concentration bottom Wall','Max Allowed Saltconcentration, 150[mg/l]')

%% Mixingzone
% output is Depth, sigma, alpha of the mixingfront
% Salitinefront;

%% Report
fprintf('Recovery starts at %f [min].\n',ConcentrationWell(1,o));
fprintf('Total amount of freshwater injected is %f[cm3].\n',max(QInjected));
fprintf('Total amount of freshwater recovered is %f[cm3], before is max salt concentration of %d[mg/l] is being reached.\n',(-1*QRetreived(1,MaxUsableCol(1))),MaxAllowedC);
fprintf('This gives a recovery efficiency of %f[percentage].\n',RE(1,MaxUsableCol));
if SizeTooSmall,
    fprintf('Size Zoetwatervat too small, leakage underneath the wall. Starting at %d[min].\n',SizeTooSmallTime);
else
    fprintf('Het zoetwatervat is large enough.\n')
end
fprintf('Injection should have been stopped at %d [min].\n',ConcentrationWell(1,MaxUsableCol));

if (ConcentrationWell(end,end)<=MaxAllowedC),
    fprintf('there is still freshwater left in vessel');
end

%% Save as file

ConcentrationSaved=XS(C(MaxUsableCol).values);
save ConcSaved ConcentrationSaved;    
fprintf('Concentration at stop point saved in "ConcSaved.mat".\n')

load results

if SizeTooSmall,
    SizeTooSmallSaved(n,1)=SizeTooSmallTime;
else
    SizeTooSmallSaved(n,1)=0;
end
ConcentrationWellSaved(n,1)=ConcentrationWell(1,m);  % After t=x[min], start retreiving
QInjectedSaved(n,1)=QInjected(m);  % Total injected
QRetreivedSaved(n,1)=(-1*QRetreived(1,MaxUsableCol(1)));  % Total retreived
RESaved(n,1)=RE(1,MaxUsableCol);  % Max recovery, untill conc exceedes maxconc
MaxQuantitySaved(n,1)=MaxQuantity;  % Max storeable freshwater

save results ConcentrationWellSaved QInjectedSaved QRetreivedSaved RESaved SizeTooSmallSaved MaxQuantity MaxQuantitySaved -append;